      Program main
      implicit none
      Integer,Parameter :: NATOM=2  !atom number
      Integer,Parameter :: NORB=10  !orbital number(AO)
      Integer,Parameter :: NMO=10   !orbital number(MO)
      Integer,Parameter :: Num_ele=14 !ele number
 
      integer i,j,k,l,ii,jj,guess
      Real*8  GEOM(NATOM,3)
      Integer ATOMCHG(Natom),iERRO
      Real*8 Hcore(NORB,NORB)       !One-electron integral 
      Real*8 S(NORB,NORB)           !Overlap matrix
      Real*8 S_DIAG(NORB,NORB)      !S^(-0.5)
      real*8 U_matrix(NORB,NORB)    !Equal to Overlap matrix
      Real*8 U(NORB,NORB,NORB,NORB) ! atom Two-electron integral
      character     :: baselable*30
      !*****************add-HF**********************
      integer :: kbasst(NATOM),kbased(NATOM),kcenter(NATOM)
      character*25 file_name
      real*8 X_transfer(NORB,NORB),X_adjoint(NORB,NORB)
      real*8 P_old(NORB,NORB) ,P_new(NORB,NORB),G_ij(NORB,NORB)
      real*8 Fock(NORB,NORB),Fock_x(NORB,NORB)
      real*8 C_prime(NORB,NORB),C_new(NORB,NORB),Cof(NORB,NORB)
      real*8 E_nu,E_ele,E2,delta
      real*8 P_den(NORB,NORB),orb_energy(NORB),grad_rhf(3*NATOM)
      real*8 Fock_xx(NORB,NORB),ING(NORB,NORB)
      real*8 J_martix(NORB,NORB),K_martix(NORB,NORB)
      integer Zn(NATOM)        !the atom charges
      logical logical_e
      integer count_iter
      !*****************Lowdin population*******
      !*****************damping and shift*******
      integer :: inv_damp
      real*8  :: damp
      integer :: inv_shift
      real*8  :: delta_e
      !*************UHF*************************
      integer :: UHF1 !!!Whether to call UHF
      integer :: n_el_a,n_el_b
      real*8  :: F_a(NORB,NORB),F_b(NORB,NORB)
      real*8  :: P_a(NORB,NORB),P_b(NORB,NORB)
      real*8  :: E_UHF,E_UHF_ele
      !******************************************
      !******************************************
      integer :: ROHF1 !!!Whether to call ROHF
      !*********************MP2******************
      real*8 Hcore_mo(NORB,NORB),U_mo(NORB,NORB,NORB,NORB)
      integer m,n,o,p,q,r,ss
      real*8  E_mp2,E_mp3
      integer Num
      Real*8 tmp1(NORB,NORB,NORB,NMO)
      Real*8 tmp2(NORB,NORB,NMO,NMO)
      Real*8 tmp3(NORB,NMO,NMO,NMO)
      Real*8 tmp4(NMO,NMO,NMO,NMO)
      real*8 t1,t2,HF_t
      !***************DIIS********************
      Integer,Parameter :: ndiis=7
      integer inv_DIIS,kq,kp,infodiis
      integer ipivdiis(ndiis+1)
      real*8 sumdiis
      real*8 fockdiis(NORB,NORB,ndiis),errmat(NORB,NORB,ndiis)
      real*8 cofdiis(ndiis+1),Bdiis(ndiis+1,ndiis+1)
      integer ierr,lwork
      real*8,allocatable :: work(:)
      real*8,allocatable :: ipiv(:)
      real*8 FPS(NORB,NORB),SPF(NORB,NORB)
      !***************End***********************
      !************** CCSD**********************
      integer :: inv_CCSD
      real*8 ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
      real*8 E_cc
      !*****************************************
      integer :: spin_multi
      !*****************************************
      call cpu_time(t1)
      !****************
      !guess=3 !!Set Fock martix initial guess 
      UHF1=0      !!=1:call UHF
      ROHF1=0     !!=1:call ROHF
      spin_multi=1
      inv_DIIS=1  !!=1:invoke DIIS
      inv_damp=0  !!=1:invoke Damp
      damp=0.2d0  !!set damp value
      inv_shift=0 !!=1:invoke level shift
      delta_e=0.2 !! set delta energy  
      inv_CCSD=1  !!=1: call CCSD
      !****************


      open(21,file='output.log')
 
      file_name='GEOM1.xyz'
      baselable='sto-3g'
      call logo(21)
      write(21,*)"Molecule and Coordinates"
      CALL RD_GEOM(NATOM,GEOM,ATOMCHG,iERRO,file_name)
      CALL INTXC_ALLOCATE(NATOM)
      CALL BASIS_INT(kbasst,kbased,kcenter,NATOM,NORB,GEOM,ATOMCHG,S,Hcore,U,baselable)
      CALL INTXC_DEALLOCATE()
      Zn=ATOMCHG
      write(21,*)'Basic set:',baselable
      write(21,*)'==========================================='
      write(21,*)'overlap matrix'
      do i=1,NORB
                write(21,*)S(i,:)
      enddo
      write(21,*)
     !write(21,*)'***************'
     !  open(26,file='U.txt')
     !do i=1,NORB
     !   do j=1,NORB
     !      do k=1,NORB
     !         do l=1,NORB
     !          WRITE(26,*)i,j,k,l,U(i,j,k,l)
     !         enddo
     !      enddo
     !   enddo
     !enddo  
     !do i=1,NORB
     !          WRITE(27,*)Hcore(i,:)
     !enddo
     !close(26)
        
      U_matrix=S

      CALL DIAG(U_matrix,S_DIAG,NORB) 

      do i=1,NORB
               S_DIAG(i,i)=S_DIAG(i,i)**(-0.5d0)
      enddo
       
      !orthogonalization transformation matrix(eigenvectors*eigenvalues)
      X_transfer = matmul(U_matrix,S_DIAG)
      
      !X_transfer transpose get X^+
      X_adjoint =transpose(X_transfer)
      call E_nuclear(NATOM,GEOM,Zn,E_nu)
      write(21,*)

      if(UHF1==1)then
        call UHF(NORB,F_a,F_b,P_a,P_b,n_el_a,n_el_b,Num_ele,spin_multi,S,Hcore,X_transfer,U,E_nu)
        !call UHFM(Hcore,U,NORB,S,E_nu)
      endif
        

      if(ROHF1==1)then
        call ROHF(NORB,X_transfer,F_a,F_b,P_a,P_b,S,Hcore,U,n_el_a,n_el_b,E_nu)
      endif

     !E_UHF=E_UHF_ele+E_nu
     !write(21,*)'Nuclear energy:      ',E_nu
     !write(21,*)'UHF Energy:  ',E_UHF
     !write(21,*)'***************************'
     !write(21,*) 

     !check: ING=x^+*S*x=I
     !ING=matmul(matmul(X_adjoint,S),X_transfer)
     !open(15,file='XSX.txt')
     !write(15,*)'this is transformed basis function orthonormal'
     !do i=1,NORB
     !   write(15,*)ING(i,:)
     !enddo
      guess=1 !!!Set Fock martix initial guess
      call  initial_guess(guess,NORB,X_adjoint,X_transfer,Hcore,S,P_old,Num_ele,E2)
      count_iter=0
      !P_old=0.0d0 
      logical_e=.true.
      !E2=0.0d0
      Fock_xx(:,:)=0.0d0 
      !========================
      !  SCF CYCLE
      !========================
      !do i = 1,1
      do while(logical_e)
        count_iter=count_iter+1

        !call G_matrix(P_old,U,NORB,G_ij)  
        call J_and_K(P_old,U,NORB,J_martix,K_martix,G_ij)
        !write(*,*)G_ij
        !stop
        Fock = Hcore + G_ij
        !!!Direct Inversion in the Iterative Subspace(DIIS)
        if(inv_DIIS==1)then !invoke DIIS
               ! write(21,*)'invoke DIIS'
                fockdiis(:,:,mod(count_iter,ndiis)+1)=Fock
            !   do i=1,ndiis
            !           write(100,*)Fock(i,:)
            !   enddo
            !   write(100,*)'===='
            !   do i=1,ndiis
            !           write(100,*)Fockdiis(i,:,2)
            !   enddo
               ! write(*,*)Fock(3,1),fockdiis(3,1,1) 
               ! stop
                !error matrix:ei=FPS-SPF
                !FPS
                FPS=matmul(matmul(Fock,P_old),S)
                !SPF
                SPF=matmul(matmul(S,P_old),Fock)
                !error_matrix=FPS-SPF
                errmat(:,:,mod(count_iter,ndiis)+1)=FPS-SPF
                !construct B matrix
                if(count_iter>=ndiis)then
                 !do i=1,ndiis
                 !      write(99,*)fockdiis(:,:,i)
                 !enddo
                  Bdiis=0.0d0
                  do kp=1,ndiis
                     do kq=1,kp
                        Bdiis(kp,kq)=sum(errmat(:,:,kp)*errmat(:,:,kq))
                        !if(kp/=kq)Bdiis(kq,kp)=Bdiis(kp,kq)
                     enddo
                  enddo
               !do i=1,ndiis+1
               !        write(99,*)Bdiis(i,:)
               !enddo
               !write(99,*)'======'
                !stop
                Bdiis(ndiis+1,:)=-1.0d0
                Bdiis(ndiis+1,ndiis+1)=0.0d0
                cofdiis=0.0d0
                cofdiis(ndiis+1)=-1.0d0
                !solve New matrix of coefficients
                lwork=-1
                allocate(ipiv(size(cofdiis)))
                do ii = 1,2
                        !lwork=-1
                        allocate(work(abs(lwork)))
                        call dsysv('L', ndiis+1, 1, Bdiis, ndiis+1, ipiv, cofdiis,ndiis+1, work, lwork, ierr)
                        lwork = nint(work(1))
                        deallocate(work)
                end do
                if(ierr/=0)stop 'Solve cofdiis failed!!!'
                deallocate(ipiv)
                !write(*,*)cofdiis
                !stop
              !call dgesv(ndiis+1,2,Bdiis,ndiis+1,ipivdiis,Cofdiis,ndiis+1,infodiis)
              !if(infodiis/=0)stop 'Solve cofdiis failed!!!'
                !set new fock matrix
                Fock=0.0d0
                do i=1,ndiis
                   Fock=Fock+cofdiis(i)*fockdiis(:,:,i)
                enddo
                endif
        endif
        if(inv_shift==1)then !!no test
                call level_shift(delta_e,NORB,S,P_old,Fock)
        endif
        Fock_x = matmul(X_adjoint,matmul(Fock,X_transfer))
        
        Fock_xx=Fock_x
        !now,Fock_prime output is an eigenvector fock_prime=C'
        call DIAG(Fock_xx,C_prime,NORB)

        !C=X*C' here,X:X_transfer C':eigenvector
        C_new=matmul(X_transfer,Fock_xx)

        !Update the P_density matrix by new coefficient(C_new)
        call P_n_den(C_new,P_new,NORB,Num_ele)
        if(inv_damp==1)then
                call damping(damp,NORB,P_new,P_old)
        endif
        !Calculated energy
        call calc_E0(Hcore,Fock,P_new,NORB,E_ele)
        !write(*,*) E_ele 
        !P_old = P_new
        call E_nuclear(NATOM,GEOM,Zn,E_nu)
        !write(*,*)Zn
        call delta_P(NORB,P_old,P_new,delta)
        if((delta < 1.0E-4).and.(abs(E_ele-E2) < 1.0E-6))then
        !if((abs(E_ele-E2) < 10.0E-6))then
               logical_e=.false.
               write(21,*)"==========================================="
               write(21,*)'Density matrix'
               do i=1,NORB
                  write(21,*)P_new(i,:)
               enddo
               write(21,*)
               write(21,*)'==========================================='
               write(21,*)'Congratulations, SCF has converged!'
               write(21,*)'Number of iterations:',count_iter
               write(21,*)'Electronic energy:   ',E_ele
               write(21,*)'Nuclear energy:      ',E_nu
               write(21,*)'Total energy(HF):    ',E_ele+E_nu
               call DIAG_1(Fock_x,orb_energy,NORB)
              ! write(22,*)'G_ij======='
             !do i=1,NORB
             ! write(22,*)P_new(i,:)
             !enddo 
               write(21,*)"Orbital Energy:      "
              do i=1,NORB
               write(21,'(2X,A3,I2,A1,3X,F16.12)')"orb",i,":",orb_energy(i)
              enddo
        endif
        P_old = P_new
        E2 = E_ele
        !write(*,*)'Number of iterations:',count_iter
        !write(*,*)E0_rhf
      enddo
      call popany(NATOM,NORB,C_new,S,Num_ele,kbasst,kbased,kcenter,ATOMCHG,P_new)
      call cpu_time(t2)
      HF_t = t2-t1
      write(21,*)'HF_TIME:',HF_t
      write(21,*)'==========================================='
      write(21,*)
!!!
        write(21,*)
        write(21,*)'****************************'
        !call geomopt_car(NATOM,Num_ele,NORB,GEOM,ATOMCHG,baselable)
        call Numgrad(NATOM,Num_ele,NORB,GEOM,ATOMCHG,baselable,grad_rhf)
        write(21,*)'****************************'
        write(21,*)
!!!
 
      
  
      !======================BEGIN OF MP2===========================      
      write(21,*)'*** Start MP2 procedure ***'
      Hcore_mo = 0.0d0
      !Calculate the one-electron integral of the molecular orbita
      do i=1,NORB  !a
        do j=1,NORB  !b
          do k=1,NORB  !u
            do l=1,NORB  !v
              Hcore_mo(i,j) = Hcore_mo(i,j)+Hcore(k,l)*C_new(k,i)*C_new(l,j)
            end do
          end do
        end do
      end do
      !Accuracy of one-electron integration of molecular orbitals
      !compared with BDF
      
     !do i=1,NORB
     !   write(23,*)Hcore_mo(i,:)
     !enddo
     !write(23,*)'*************'
      !========================================================

      !========================================================
      !Calculate the one-electron integral of the molecular orbita
      U_mo = 0.0d0
      !该方法计算量N**8 
      do i=1,NORB
        do j=1,NORB
          do k=1,NORB
            do l=1,NORB
      !-----m,n,o,p equal to the atomic orbital integral-----------
              do m=1,NORB
                do n=1,NORB
                  do o=1,NORB
                    do p=1,NORB
                    U_mo(i,j,k,l)= U_mo(i,j,k,l)+U(m,n,o,p)*C_new(m,i)*C_new(n,j)*C_new(o,k)*C_new(p,l)
                    end do
                  end do
                end do
              end do
            end do
          end do
       end do
      end do
      
      !=====================================================
      !Accuracy of TWO-electron integration of molecular orbitals
      !compared with BDF
      !=====================================================
     !do i=1,NORB
     !   do j=1,NORB
     !          do k=1,NORB
     !                  do l=1,NORB
     !                   write(23,*)U_mo(i,j,k,l),'   ',i,j,k,l
     !                  enddo
     !          enddo
     !  enddo
     !enddo
     !======================================================
       tmp1=0.0d0
      do i=1,NORB
        do j=1,NORB
           do k=1,NORB
              do ss=1,NMO
                 do l=1,NORB
                  !C_new(l,ss)*(ij|kl) => (ij|ks)
                  tmp1(i,j,k,ss)=tmp1(i,j,k,ss) + C_new(l,ss)*U(i,j,k,l)
                 enddo
              enddo
           enddo
        enddo
      enddo
      tmp2=0.0d0
      do i=1,NORB
        do j=1,NORB
           do r=1,NMO
              do ss=1,NMO
                 do k=1,NORB
                  !C_new(l,ss)*(ij|ks) => (ij|rs)
                  tmp2(i,j,r,ss)=tmp2(i,j,r,ss) + C_new(k,r)*tmp1(i,j,k,ss)
                 enddo
              enddo
           enddo
        enddo
      enddo
      tmp3=0.0d0
      do i=1,NORB
        do q=1,NMO 
           do r=1,NMO
              do ss=1,NMO
                 do j=1,NORB
                  !C_new(l,ss)*(ij|rs) => (iq|rs)
                  tmp3(i,q,r,ss)=tmp3(i,q,r,ss) + C_new(j,q)*tmp2(i,j,r,ss)
                 enddo
              enddo
           enddo
        enddo
      enddo
      tmp4=0.0d0
      do p=1,NMO
        do q=1,NMO 
           do r=1,NMO
              do ss=1,NMO
                 do i=1,NORB
                  !C_new(l,ss)*(iq|rs) => (pq|rs)
                  tmp4(p,q,r,ss)=tmp4(p,q,r,ss) + C_new(i,p)*tmp3(i,q,r,ss)
                 enddo
              enddo
           enddo
        enddo
      enddo
     !========================================================
      U_mo=tmp4
                     
     !======================================================
      E_mp2 = 0.0d0
      Num = Num_ele/2
      
      do i=1,Num    !a
        do j=1+Num,NORB !r
          do k=1,Num       !b
            do l=1+Num,NORB   !s
            E_mp2 = E_mp2 + U_mo(i,j,k,l)*(2*U_mo(i,j,k,l)-U_mo(i,l,k,j))/(orb_energy(i)+orb_energy(k)-orb_energy(j)-orb_energy(l))
       
            enddo
          enddo
        enddo
      enddo

      write(21,'(35X,40A)')'*******************'
      write(21,'(35X,40A)')'*       MP2       *'
      write(21,'(35X,40A)')'*******************'

      write(21,*)'E_mp2:            ', E_mp2
      write(21,*)'Total energy(mp2):',(E_nu+E_ele)+E_mp2
      write(21,*)'*** The MP2 procedure Done **'
      write(21,*)'==========================================='
      write(21,*)
      !!!=============BEGINING OF MP3=======================
      write(21,*)'*** Start  MP3 procedure ***'
      call MP3(Num_ele,NORB,orb_energy,U_mo,E_mp3)
      write(21,'(35X,40A)')'*******************'
      write(21,'(35X,40A)')'*       MP3       *'
      write(21,'(35X,40A)')'*******************'

      write(21,*)'E_mp2:            ', E_mp3
      write(21,*)'Total energy(mp3):',(E_nu+E_ele)+E_mp2+E_mp3
      write(21,*)'*** The MP3 procedure Done **'
      write(21,*)'==========================================='
      write(21,*)
      !!!END MP3
        
      !===============Begining of CCSD=======================
      E_mp2=0.0d0
      !write(1000,*)U_mo(1,1,1,1)
      call MO2SOERI(NORB,U_mo,ERISPIN)  !MO to spin-orbitalERI(chem)<pq||rs>
      !call MO2SOERI1(NORB,U_mo,ERISPIN)!! phy 
      !!!===============================================================
      !!!Spin orbital integral is used to calculate E_MP2
      call mp2_energy(Num_ele,NORB,orb_energy,ERISPIN,E_mp2) 
      write(21,*)"Spin orbit integral is used to calculate E_MP2:",E_mp2 
      !!!================================================================
      if(inv_CCSD==1)then
      call CC(E_ele,E_nu,Num_ele,NORB,orb_energy,ERISPIN,E_cc)
      endif 
      !!!END CCSD(T)
      !!!====================================================
      !!!Test  Slater-Condon rules 
      !!!The ground state energy of the molecule is calculated using the
      !spin orbit integra<ab||ab>
      call e1(Num_ele,Hcore_mo,ERISPIN,NORB) 
      !!! Test correctly
      !!!===============Begining of CISD=======================
      close(21)
      ENDPROGRAM

      SUBROUTINE RD_GEOM(NATOM,GEOM,ATOMCHG,iERRO,file_name)
      IMPLICIT NONE
      INTEGER NATOM,iERRO,I
      INTEGER ATOMCHG(NATOM)
      character*25 file_name
      REAL*8  GEOM(NATOM,3)
      CHARACTER*10 CTYPE(NATOM)
      GEOM=0.0D0
      ATOMCHG=0
      OPEN(23,FILE=file_name)
      !write(*,*)file_name
      DO I=1,NATOM
        READ(23,*)CTYPE(I),GEOM(I,1),GEOM(I,2),GEOM(I,3)
        write(21,*)CTYPE(I),GEOM(I,1),GEOM(I,2),GEOM(I,3)
      ENDDO
      CLOSE(23)
      DO I=1,NATOM
        IF(CTYPE(I)=='H'  .or. CTYPE(I)=='h')  ATOMCHG(I)=1
        IF(CTYPE(I)=='He' .or. CTYPE(I)=='he') ATOMCHG(I)=2
        IF(CTYPE(I)=='Li' .or. CTYPE(I)=='li') ATOMCHG(I)=3
        IF(CTYPE(I)=='Be' .or. CTYPE(I)=='be') ATOMCHG(I)=4
        IF(CTYPE(I)=='B'  .or. CTYPE(I)=='b')  ATOMCHG(I)=5
        IF(CTYPE(I)=='C'  .or. CTYPE(I)=='c')  ATOMCHG(I)=6
        IF(CTYPE(I)=='N'  .or. CTYPE(I)=='n')  ATOMCHG(I)=7
        IF(CTYPE(I)=='O'  .or. CTYPE(I)=='o')  ATOMCHG(I)=8
        IF(CTYPE(I)=='F'  .or. CTYPE(I)=='f')  ATOMCHG(I)=9
        IF(CTYPE(I)=='Ne' .or. CTYPE(I)=='ne') ATOMCHG(I)=10
      ENDDO
      END SUBROUTINE RD_GEOM

      SUBROUTINE DIAG(A,S_DIAG,NORB) 
       real*8  A(NORB,NORB),S_DIAG(NORB,NORB)
       real*8  W(NORB),V_m(NORB,NORB)
       integer, parameter :: lwmax=10000
       real*8, dimension(lwmax):: work
       integer i,j,N,LDA,lwork,info

       
       LDA = NORB
       N = NORB
       lwork = -1
       !allocate (work(lwork))
       call dsyev('V', 'U', N, A, LDA, W, work, lwork,info)
       lwork =min(lwmax,int(work(1)))*2*N
       !allocate (work(lwork))
       !call dgeev('V','V',N,A,LDA,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
       call dsyev('V', 'U', N, A, LDA, W, work, lwork,info)
      
       if (INFO .NE. 0) THEN
       write(*,*) "ERROR: IMPOSSIBLE TO SOLVE THE EIGENVALUE PROBLEM!" 
       endif

       S_DIAG = 0
       do i=1,NORB
                S_DIAG(i,i)=W(i)
       enddo
      !open(23,file='S_diag')
      ! write(23,*)'Diagonalize the over matrix'
      !do i=1,NORB
      !         write(23,*)S_DIAG(i,:)
      !         write(24,*)A(i,:)
      !enddo
       ENDSUBROUTINE
        
       SUBROUTINE DIAG_1(A,W,NORB)
       real*8  A(NORB,NORB)
       real*8  W(NORB),V_m(NORB,NORB)
       integer, parameter :: lwmax=10000
       real*8, dimension(lwmax):: work
       integer i,j,N,LDA,lwork,info
       LDA = NORB
       N = NORB
       lwork = -1
       !allocate (work(lwork))
       call dsyev('V', 'U', N, A, LDA, W, work, lwork,info)
       lwork =min(lwmax,int(work(1)))*2*N
       !allocate (work(lwork))
       !call
       !dgeev('V','V',N,A,LDA,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
       call dsyev('V', 'U', N, A, LDA, W, work, lwork,info)

       if (INFO .NE. 0) THEN
       write(*,*) "ERROR: IMPOSSIBLE TO SOLVE THE EIGENVALUE PROBLEM!"
       endif
       
      !open(23,file='S_diag')
      ! write(23,*)'Diagonalize the over matrix'
      !do i=1,NORB
      !         write(23,*)S_DIAG(i,:)
      !         write(24,*)A(i,:)
      !enddo
       ENDSUBROUTINE

      
      !Inverse matrix  
      subroutine INV_MATRIX(m,inva,Nr)
        implicit none
        integer :: i,j
        integer :: info,Nr
        integer :: ipiv(Nr)
        real*8  :: work(Nr)
        real*8 :: m(Nr,Nr),inva(Nr,Nr),n(Nr,Nr)
        !write(*,*)Nr
        !do i =1,Nr           
        !        write(13,*)m(i,:)
        !enddo
        call dgetrf(Nr,Nr,inva,Nr,ipiv,info)
        call dgetri(Nr,inva,Nr,ipiv,work,Nr,info)
        n=matmul(m,inva)
        !do i =1,Nr
        !       write(13,*)n(i,:)
        !enddo
        end subroutine
        
        subroutine J_and_K(P_density,U,NORB,J_martix,K_martix,G_ij)
        implicit none
        !============in===========
        integer i,j,k,l
        integer NORB
        real*8 P_density(NORB,NORB)
        real*8 U(NORB,NORB,NORB,NORB)
        !=============out==========
        real*8 J_martix(NORB,NORB),K_martix(NORB,NORB)
        real*8 G_ij(NORB,NORB)

        J_martix(:,:) = 0.0d0
        K_martix(:,:) = 0.0d0
        do i =1,NORB
           do j=1,NORB
              do k=1,NORB
                 do l=1,NORB
                     J_martix(i,j) = J_martix(i,j) + P_density(k,l)*U(i,j,l,k)
                     K_martix(i,j) = K_martix(i,j) + 0.5d0*P_density(k,l)*U(i,k,l,j) 
                 enddo
              enddo
           enddo
        enddo
        do i=1,NORB
           do j=1,NORB
           G_ij(i,j) = J_martix(i,j) - K_martix(i,j)
           enddo
        enddo
       endsubroutine 
        
        subroutine G_matrix(P_density,U,NORB,G_ij)
        integer i,j,k,l,NORB
        real*8 U(NORB,NORB,NORB,NORB),G_ij(NORB,NORB)
        real*8 P_density(NORB,NORB)
        real*8 temp
        real*8 doublej,doublek
        G_ij(:,:)=0.0d0
        
        do i=1,NORB
           do j=1,NORB
              do k=1,NORB
                 do l=1,NORB
                        doublej=U(i,j,l,k)
                        doublek=0.5d0*U(i,k,l,j)            
                        G_ij(i,j)=G_ij(i,j) + P_density(k,l)*(doublej-doublek)
                        
                 enddo
              enddo
     
           enddo
        enddo
        endsubroutine   
        

        subroutine P_n_den(C_new,P_new,NORB,Num_ele)
        !input
        integer NORB
        real*8 C_new(NORB,NORB)       !New coefficient matrix
        integer i,j,k,x
        integer Num_ele
        
        !output
        real*8 P_new(NORB,NORB)       !New density matrix
        
        x=Num_ele/2
        !write(*,*)x
        P_new=0.0d0
        
        do i=1,NORB
                do j=1,NORB
                   do k=1,x
                      P_new(i,j)=P_new(i,j)+2.0d0*C_new(i,k)*C_new(j,k)
                   enddo
                enddo
        enddo
        endsubroutine P_n_den
        
        subroutine calc_E0(Hcore,Fock,P_new,NORB,E_ele)
        !input
        real*8 Hcore(NORB,NORB),Fock(NORB,NORB)
        real*8 P_new(NORB,NORB)
        integer NORB
        
        !output
        real*8 E_ele
        integer i,j
        
        E_ele=0.0d0
        do i=1,NORB
                do j=1,NORB
                    E_ele =E_ele + 0.5d0*P_new(j,i)*(Hcore(i,j)+Fock(i,j))
                enddo
        enddo
        !write(*,*)E0_rhf
        endsubroutine calc_E0
        
        subroutine E_nuclear(Nn,Rn,Zn,E_nu) 
        implicit none

        ! input
        integer Nn            ! Totalnumber of nuclei
        integer Zn(Nn)        ! Nuclear charges
        real*8  Rn(Nn,3)      ! Nuclear positions
        real*8  a             ! bohr     
        integer :: i,j        ! Loop variables
        ! output
        real*8  E_nu          ! Nuclear energy

        E_nu = 0.0D0
        a = 0.5291772083 !bohr= ai*a
        do i = 1, Nn
            do j = i + 1, Nn
                E_nu = E_nu + Zn(i)*Zn(j)/NORM2(Rn(i,1:3)/a-Rn(j,1:3)/a)
            enddo
        enddo
            
        endsubroutine E_nuclear      


        subroutine delta_P(NORB,P_old,P_new,delta) 
        implicit none

        ! input
        integer  NORB          
        real*8  P_old(NORB,NORB)    
        real*8  P_new(NORB,NORB)    
        integer :: i,j

       ! output
       ! Sum ofmatrix elements square differences
        real*8 :: delta      !Second convergence criterion 
        

        delta = 0.0d0
        do i = 1, NORB
            do j = 1, NORB
                delta = delta + (P_old(i,j)-P_new(i,j))**2
            enddo
        enddo
        !dsqrt:Square root
        delta = DSQRT(delta / NORB**2)

        endsubroutine delta_P

        subroutine mp2_energy(Num_ele,NORB,Eorb,ERISPIN,E_mp2tem)
        implicit none
        integer :: NORB
        integer :: Num_ele,nspin,nspinocc,nspinvir
        integer :: ii,ia,ij,ib !spin orb
        integer :: iispa,iaspa,ijspa,ibspa
        real*8  :: Eorb(NORB)
        real*8  :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        real*8,allocatable::tijab(:,:,:,:),Dijab(:,:,:,:)
        real*8  :: E_mp2tem,tijabnow
        nspin=2*NORB
        nspinocc=Num_ele
        nspinvir=nspin-nspinocc
        !write(1000,*)"3",ERISPIN(1,2,1,2)
        allocate(tijab(nspinocc,nspinocc,nspinvir,nspinvir),&
        Dijab(nspinocc,nspinocc,nspinvir,nspinvir))
        E_mp2tem=0d0;Dijab=0.0d0;tijabnow=0.0d0
        do ii=1,nspinocc
            iispa=ceiling(ii/2d0)
            do ia=nspinocc+1,nspin
                iaspa=ceiling(ia/2d0)
                !Dia(ii,ia-nspinocc)= Eorb(iispa)-Eorb(iaspa)
                do ij=1,nspinocc
                ijspa=ceiling(ij/2d0)
                    do ib=nspinocc+1,nspin
                    ibspa=ceiling(ib/2d0)
                    Dijab(ii,ij,ia-nspinocc,ib-nspinocc)=Eorb(iispa)+Eorb(ijspa)-Eorb(iaspa)-Eorb(ibspa)
                    tijabnow=ERISPIN(ii,ij,ia,ib)/Dijab(ii,ij,ia-nspinocc,ib-nspinocc)
                    !write(*,*)tijabnow
                    E_mp2tem=E_mp2tem+ERISPIN(ii,ij,ia,ib)*tijabnow
                    tijab(ii,ij,ia-nspinocc,ib-nspinocc)=tijabnow
                    end do
                end do
            end do
        end do
        E_mp2tem=0.25D0*E_mp2tem
        !write(*,*)E_mp2tem
        endsubroutine mp2_energy

        subroutine popany(NATOM,NORB,Cof,S,Nele,kbasst,kbased,kcenter,ATOMCHG,P)
        implicit none
        integer :: i,j,nocc
        integer :: NATOM,Ntem,Nele,NORB
        real*8 :: Cof(NORB,NORB),S(NORB,NORB)
        real*8 :: Den(NORB,NORB),pop(NORB,NORB),P(NORB,NORB)
        integer :: kbasst(NATOM),kcenter(NATOM),kbased(NATOM)
        integer :: ATOMCHG(NATOM)
        real*8  ::Zpop(NATOM)
        real*8  :: ATOMCHG1(NATOM)
        real*8 :: tem,zmu
        den=0d0
        nocc=nele/2
        if(nocc==0)nocc=1
        !cof0=transpose(Cof)
        do i=1,NATOM
                ATOMCHG1(i)=ATOMCHG(i)
        enddo
        !write(*,*)'000',ATOMCHG1(1)+1.1
        !!!method 1 
        Den=matmul(cof(:,1:nocc),transpose(cof(:,1:nocc)))
        Pop=matmul(Den,S)
        !!!end method
        !!!method 2
        !**Pop=matmul(P,S)
        !!!end method
        !write(*,*)Pop
        Ntem=1
        tem=0d0
        do i=1,NATOM
            if(kcenter(i)==Ntem)then
                do j=kbasst(i),kbased(i)
                    tem=tem+2*pop(j,j) !!!method1
                    !tem=tem+pop(j,j)  !!!method2
                end do
            else
                !write(*,*)tem
                Zpop(kcenter(i-1))=ATOMCHG1(kcenter(i-1))-tem
                !write(*,*)'222',ATOMCHG1(kcenter(i-1)),tem
                !Zpop(kcenter(i-1))=ATOMCHG1(kcenter(i-1))-tem
                !write(*,*)'111',Zpop(kcenter(i-1))
                ntem=kcenter(i)
                tem=0d0
                do j=kbasst(i),kbased(i)
                     tem=tem+2*pop(j,j) !!!method 1
                     !tem=tem+pop(j,j)  !!!method 2
                end do
                !write(*,*)tem
            end if
            if(i==NATOM)then
                ZPop(kcenter(i))=ATOMCHG1(kcenter(i))-tem
            end if
        end do
        write(21,*)'==========================================='
        write(21,*)'*** Population Analysis/Atomic Charges ***'
        write(21,*)'               Atomic charges',&
                        '            Mulliken charges on atom'
        do i=1,natom
            write(21,*)'ATOM',i,ATOMCHG(i)-ZPop(i),ZPop(i)
        end do
        ZMu=sum(ZPop)
        if(abs(ZMu)<1d-8)zmu=0d0
        write(21,*)'Sum of Mulliken charges =  ',zmu
        write(21,*)'==========================================='
        endsubroutine popany
        
        subroutine e1(Num_ele,Hcore,ERISPIN,NORB)  !!! Ground state energy
        implicit none
        integer :: Num_ele,i,j,ia,NORB
        real*8  :: haa,h22,e11
        real*8  :: Hcore(NORB,NORB),ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
        e11=0.0d0;haa=0.0d0;h22=0.0d0
        do i=1,Num_ele
           ia=ceiling(i/2d0)
           haa=Hcore(ia,ia)
           h22=0.0d0
           do j=i,Num_ele
                 h22=h22+ERISPIN(i,j,i,j)
           enddo
          e11=haa+h22+e11
        enddo
        !write(*,*)e11
       endsubroutine 
